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Abstract 

We formulate an optimal load control (OLC) problem in power networks where the objective 
is to minimize the aggregate cost of tracking an operating point subject to power balance over the 
network. We prove that the swing dynamics and the branch power flows, coupled with frequency-based 
load control, serve as a distributed primal-dual algorithm to solve OLC. Even though the system has 
multiple equilibrium points, we prove that it nonetheless converges to an optimal point. This result 
implies that the local frequency deviations at each bus convey exactly the right information about the 
global power imbalance for the loads to make individual decisions that turn out to be globally optimal. It 
allows a completely decentralized solution without explicit communication among the buses. Simulations 
show that the proposed OLC mechanism can resynchronize bus frequencies with significantly improved 
transient performance. 

A preliminary version has appeared in the Proceedings of the 3 rd IEEE International Conference on Smart Grid Communica- 
tions, Tainan City, Taiwan, November 2012. 
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I. Introduction 

In power systems, regulation efforts traditionally focus on the generation side. For example, 
the automatic generation control adjusts the setpoints of generators based on area frequency 
deviations and unscheduled cross-area power flows |0Q|. To track the given setpoints at a faster 
time scale and improve stability, control mechanisms, such as the excitation system, speed 
governing system and power system stabilizer, are deployed on the generation side [EJ [I3J . 
However, relying solely on generation control may not be enough. Due to limited ramping 
rate and large inertia, generators are suitable for minute-by-minute power balance, but may 
incur expensive wear-and-tear, high emissions, and low thermal efficiency when responding 
to regulation signals at intervals of seconds ifl[|5]|. Complementary to generation control, we 
consider load control as an additional mechanism that provides fast and inexpensive power system 
regulation. Indeed the feasibility and efficiency of load control has already been demonstrated 
in several electricity markets. Long Island Power Authority (LIPA) developed LIPAedge, which 
provides 24.9 MW of demand reduction and 75 MW of spinning reserve by 23,400 loads O. 
Electric Reliability Council of Texas (ERCOT) has 50% of its 2400 MW reserve provided by 
loads. Pennsylvania- New Jersey-Maryland Interconnection (PJM) opens up the regulation market 
to participation by loads [4J. While most of the existing programs focus on direct manipulation of 
loads in a centralized scheme, the alternative strategy of decentralized load control via frequency 
measurement has also been studied in the literatures. Brooks et al. suggests that loads can sense 
and respond to frequency and provide regulation within 1 second flU. Molina-Garcia et al. 
studies the aggregate response characteristics when individual loads are turned on/off as the 
frequency fluctuates [7J. Donnelly et al. develops proportional control of intelligent loads, and 
investigates the effect of distribution systems, the effect of discretized control, and the effect of 
time-delay of control actions, through the simulation of a 16-generator transmission network (8]|. 
In our previous papers fl9ll— ffTTTl we study a decentralized algorithm to minimize the cost of load 
control based on local frequency measurement and neighborhood communication. Frequency- 
based load control does not require communication to a centralized grid operator, and is thus 
suitable for large-scale decentralized deployment. 

In this paper we consider a power transmission network where, at steady state, the generator 
frequencies at different buses (or in different balancing authorities) are synchronized to the 
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same nominal value and the mechanic power is balanced with the electric power at each bus. 
Suppose a small change in generation occurs on an arbitrary subset of the buses. How should 
the controllable loads in the network be reduced (or increased) in real time in a way that (i) 
balances the generation shortfall (or surplus), (ii) resynchronizes the bus frequencies, and (iii) 
minimizes a measure of aggregate cost of participation in such a load control? We formalize 
this question as an optimal load control (OLC) problem. The basic dynamics at each generation 
bus is described by swing equations that relate the imbalance between generation and load to 
the rate of frequency change. We assume the change in generation is small and the DC load 
flow model is reasonably accurate. We develop a frequency-based load control where loads are 
controlled based on locally measured frequency deviations and their individual cost functions. 
We prove that this frequency-based load control coupled with the system dynamics serve as 
a distributed primal-dual algorithm to solve OLC. Simulation of the IEEE 68-bus test system 
shows that the proposed OLC mechanism resynchronize bus frequencies with smaller steady- 
state error, smaller overshoot and shorter settling time, compared to the case with only local 
generator control mechanisms like power system stabilizers. 

The paper is organized as follows. Section [n] describes a dynamic model of power networks 



and formulates OLC. Section III explains how the frequency-based load control and the system 
dynamics serve as a distributed primal-dual algorithm to solve OLC and its dual. Section 
IV provides the convergence analysis of the primal-dual algorithm. Section [V] illustrates the 



proposed scheme through the simulation of the IEEE 68-bus test system. Section VI concludes 



the paper. Section VII provides a detailed justification of our model of the branch power flow, 
both analytically and through simulation of the much more detailed dynamic model of power 
systems developed in [|3l [[T2l . 



II. Problem formulation 

Let K denote the set of real numbers and N denote the set of non-zero natural numbers. For 
a set M, let \M\ denote its cardinality. A variable without a subscript usually denotes a vector 
with appropriate components, e.g., to = (u}j,j G AT) G M'-^L For a, b G E, a < b, the expression 
[■] b a denotes max{min{-, b}, a}. For a matrix A, let A T denote its transpose. For a signal u(t) 
of time, let to denote its time derivative 
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A. Transmission network model 

The power transmission network is described by a graph (J\f, £) where J\f = {1, . . . , \Af\} is 
the set of buses and £ C J\[ x M is the set of transmission lines connecting the buses. We make 
the following assumptions: Q 

• The lines (i, j) E £ are lossless and characterized by their reactance Xij. 

• The voltage magnitudes \Vj\ of buses j E M are constants. 

• Reactive power injections at the buses and reactive power flows on the lines are ignored. 
We assume that (Af,£) is directed, with an arbitrary orientation, so that if E £ then 

E~ £. We use and i — >■ j interchangeably to denote a link in £, and use "z : i — > j" 
and ll k : j — >■ A;" respectively to denote the set of buses 2 that are predecessors of bus j and 
the set of buses k that are successors of bus j. We also assume without loss of generality that 
(A/", £) is connected. 

The network has two types of buses: generator buses and load buses. A generator bus not only 
has loads, but also an AC generator that converts mechanic power into electric power through a 
rotating prime mover. A load bus has only loads but no generator. We assume that the system is 
always under the three-phase balanced condition, and for a bus j E Af, its phase a voltage at time 
t is y/2\Vj\ cos(u°t + 9® + A9j(t)), where lu° is the nominal frequency, 9°- is the nominal phase 
angle, and A9j(t) the time- varying phase angle deviation. The frequency at bus j is defined as 
ujj := cu° + A9j, and we call AtOj := A9 j the frequency deviation at bus j. We assume that 
the frequency deviations Aluj are small for all the buses j E M and the differences A9i — A9j 
between phase angle deviations are small across all the links E £. We adopt a standard 
dynamic model, e.g., in [2, Section 11.4]. 

Generator buses. We assume coherency between the internal and terminal (bus) voltage phase 
angles of the generator so that these angles always differ by a constant even during the transient, 



which is discussed in more detail in Section VII-C Then the dynamics on a generator bus j is 
modeled by the swing equation 

MjAuj + = / , r , -/ , ;,L ; -/ , ;• 

'These assumptions are similar to the standard DC approximation except that we do not assume the nominal phase angle 
difference is small across each link. 
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where Mj > is the inertia constant of the generator. The term D'-Aujj with D'- > represents 
the deviation in generator power loss due to friction, from its nominal value P® oss j . Here P™' 
is the mechanic power injection to the generator, and P? is the electric power export of the 
generator, which equals the sum of loads on bus j and the net power injection from bus j to 
the rest of the network. 

In general, load power may depend on both the bus voltage magnitude (which is assumed fixed) 
and frequency. We distinguish between two types of loads, frequency-sensitive and frequency- 
insensitive. We assume frequency- sensitive (e.g., motor-type) loads increase linearly with fre- 
quency deviation and model the aggregate of these loads by trj + D" Aujj with D" > 0, where 
d° is its nominal power. We assume frequency-insensitive loads can be actively controlled and 
our goal is to design and analyze these control laws. Let dj denote the aggregate of frequency- 
insensitive loads on bus j. Finally we allow constant power load Pj that is switched on or off. 
Then the electric power P? is the sum of frequency- sensitive loads, frequency-insensitive loads, 
constant power load, and the net power injection from bus j to other buses: 

P/ := dV + D'jAuj + dj + pi+J^Pjk-Y,^ 

k:j— >k 

where Pj k is the branch power flow from bus j to bus k. 
Hence the dynamics on a generator bus j is 

MjAuj = — [DjAuj + dj — PJ 1 + P° ut — Pj n ) 

where D 3 := + D>>, Ff := Pf - P% loss - d] - Pj, and P° ut := /> , and Pj n : = 

are respectively the total branch power flows out and into bus j. Let Pj 71 ' , P?- 
denote the nominal (operating) point at which — P™' + P° ut — APj n '° = 0, and let dj(t) = 
d°j + Adj(t), P™(t) = Pp° + APp(t), Pij{t) = PP. + APy(t). Then the deviations satisfy 

MjAujj = - (DjAcoj + Adj - APf + AP° ut - APj n ) . (1) 

Figure [TJ is a schematic diagram of the generator bus model ([T]). 

Load buses. For a load bus that has no generator, there is the following algebraic relation 
between the variables introduced above: 

= DjAujj + Adj - APp + AP° ut - APj n . (2) 

2 There may be load buses with large inertia that can be modeled by swing dynamics {TJ as proposed in [13|. We will treat 
them as generator buses mathematically. 
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for all k s.t. j [ — > k 

Fig. 1. Schematic diagram of a generator bus j, where Auij is the frequency deviation; AP™ is the change in mechanic 
power; DjAcjj characterizes the effect of generator friction and frequency-sensitive loads; Ad, is the change in aggregate 
frequency-insensitive load; APij is the deviation in branch power injected from another bus i to bus j; APjk is the deviation 
in branch power delivered from bus j to another bus fc. 



Branch flows. The deviations AP^ from the nominal branch flows follow the dynamics 

AP^ = Bij {Acoi - Acoj) , (3) 

where 

B l3 := sS^cos^ -^) (4) 

Xij 

is a constant related to nominal bus voltages and the line reactance. The same model was given in 
the literature, e.g., 01 EO, based on the quasi-steady-state assumptions. However, we derive this 
model by solving the differential equation that characterizes the line induction to obtain three- 
phase instantaneous power flow, without explicitly using the quasi-steady-state assumptions. The 



detailed derivation is given in Section VII-A 



Dynamic network model. We denote the set of generator buses by Q, the set of load buses 
by C, and use \Q\ and \C\ to denote the number of generator buses and load buses, respectively. 
Without loss of generality, let Q = {1,...,\Q\} and C = {\Q\ + 1, \M\}. In summary, the 
dynamic model of the transmission network is specified by ([T])-([3]). To simplify notation, we 
drop the A from the variables denoting deviations and write ([T])-([3]) as 



ujj = -^(DjUjj + dj-Pp + P^-PJ 1 ) for allied (5) 

)OUt 



= DjUjj + dj - Pp + P° ut - Pf for all j e C (6) 

Pij = B^ (ui - ujj) for all e £ , (7) 

where B^ are given by (Q. Hence for the rest of this paper all variables represent deviations 
from their nominal values. We will refer to the term DjUij as the deviation in the (aggregate) 
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frequency-sensitive load even though it also includes the deviation in generator power loss due 
to friction. We will refer to the term Pj™ as the change in generation even though it may also 
include changes in constant power loads. 

A steady state of the dynamic system described by ©-([7]) is defined as a state in which all 
power deviations and frequency deviations are constant over time, cbj = for j G Q, and Py = 
for e £. 

Discussion. The model given by ©-([7]) captures the power system behavior at the timescale 
of seconds. In this paper we will only consider a step change in generation (constant deviations 
Pp), which implies that the proposed model does not include the action of turbine-governor 
(at a similar or slower timescale than the proposed model) that changes the mechanic power 
injection in response to frequency deviation to rebalance power. Nor does it include any secondary 
frequency control mechanism such as automatic generation control that operates at a much slower 
timescale to restore the nominal frequency. This model therefore explores the feasibility of the 
active control of frequency-insensitive loads at a fast timescale as a supplement to the turbine- 
governor mechanism to resynchronize frequency and rebalance power. Our results in Sections 



III and M suggest it is feasible. 



The proposed model implicitly assumes that the time for individual frequencies Uj to resyn- 
chronize (converge to a common system frequency, which may be time varying) after a power 
imbalance can be similar to the time for them to converge to their equilibrium value (which may 
be different from the nominal value). Whether this assumption holds depends on the electrical 
distances 03) between different buses. For buses that are close, they resynchronize almost 
instantly to a common frequency which then converges more slowly to the equilibrium value. For 
buses that are far away, their resynchronization times are similar to convergence times. These 
observations are shown by the simulation of a more realistic model; see Figure [9] in Section 
IVIFCl 



B. Optimal load control 

Suppose a step change P m = (Pj 71 , j G AT) in generation is injected to the set J\f of buses. How 
should the frequency-insensitive loads d = (dj,j G AT) in the network be reduced (or increased) 
in real-time in a way that (i) balances the generation shortfall (or surplus), (ii) resynchronizes 
the bus frequencies, and (iii) minimizes a measure of aggregate disutility of participation in such 
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a load control? We now formalize these questions as an optimal load control (OLC) problem. 

A change P' m in generation causes a nonzero frequency deviation Uj at bus j. This frequency 
deviation incurs a cost to the aggregate frequency- sensitive load dj := DjUj and suppose this 
cost is quadratic in frequency deviations, i.e., i^yd?- = \DjUj 2 - at bus j. Suppose the aggregate 
frequency-insensitive load at bus j is to be changed by an amount dj and this change incurs a 
cost (disutility) of Cj(dj). We assume — oo < dj < dj < dj < oo. Our goal is to minimize the 
total cost over d and d while balancing generation and load across the network, written as 
OLC: 



subject to V (dj + dj) = VP™. (9) 



mm 

d<d<d,d 



Remark 1. Note that (|9]) does not require the balance of generation and load at each individual 
bus, but only balance across the entire network. This constraint is less restrictive and offers 
more opportunity to minimize costs. Additional constraints can be imposed if it is desirable that 
certain buses balance their own supply and demand, e.g., for economic or regulatory reasons. 

We assume the following condition throughout the paper: 

Condition 1. The OLC problem is feasible, and the cost functions Cj are strictly convex and 
twice continuously differentiable on [d -,^-]. 

Families of cost functions that satisfy Conditions [T] have been used in the literature [fT5l[[T6l . 
The choice of cost functions can be based on physical characteristics of the loads and user 
comfort levels. Examples of particular cost functions can be found in ifPTI for air conditioners 
and in [fT8l for plug-in electric vehicles. 



III. Load control and system dynamics as primal-dual algorithm 



We present the main results in this section, and prove them in Section IV 



A. Main results 

The objective function of the dual problem of OLC is 



E $ ^ z/ ) := E min - ( c ^ " vd i + ^7r4 _ + vP T ) ' 

• kt ■ »/■ d-<di<dj,di \ Z,Uj I 
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where can be written as 

*,-(!/) := Ci (d i (^))-^(v)-^ 2 + vP/ 1 (10) 

with 

r / , . nij 

(ID 



This objective function has a scalar variable v and is not separable across buses j G A/". Its 
direct solution hence requires coordination across buses. We propose the following distributed 
version of the dual problem over the vector v — G AT), where each bus j optimizes over 
its own variable Uj which are constrained to be equal at optimality. 
DOLC: 



max $(v) := ^ ^ 
subject to Vi = Uj for all G S. 
We have the following two lemmas regarding DOLC and its relation with OLC. 
Lemma 1. The objective function $ of DOLC is strictly concave over M)^ . 
Lemma 2. 1) DOLC has a unique optimal point v* with v* = v* = v* for all i, j G J\f. [j 



2) OLC has a unique optimal point (d*, d*) where d* = dj(u*) is given by ( fTT) and d* = D^v* 
for all j G N. 



Lemma [T] and Lemma |2] are respectively proved in Sections VIII-A and VIII-B 



Instead of solving OLC directly, Lemma [2] suggests solving its dual DOLC and recovering 
the unique optimal point (d*, d*) of the primal problem OLC from the unique dual optimal v* . 
To derive a distributed solution for DOLC, consider its Lagrangian 

H^^) ■= ^2®j(yj)~ ^2 Kijivi-Pj), (12) 
jeAf (ij)ee 

where v G R w is the (vector) variable for DOLC and 7T G is the associated dual variable 
for the dual of DOLC. Hence 7Ty, for all (i, j) G £ , measure the cost of not synchronizing the 

3 For simplicity, we abuse the notation and use u* to denote both the vector [y* , j € Af) and the common value of its 
components. Its meaning should be clear given the context. 
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variables v { and Uj across buses i and j. Using pO])-([T2]), a partial primal-dual algorithm for 
DOLC takes the form 



dL 



v 



- 7;^>, tt) = - 7j + D jVj - Pf + Trf - tt" 1 ) for j G £ (13) 
8T 

= — (z/,7r) = -(rf J (^) + J D^ J -Pf + 7rf -7T-) for j E C (14) 

7r»i = — ( z/ ' 7r ) = fm(i,j)e£, (15) 

where 7, > 0, £y > are stepsizes and 7r° ut := Y.k-.j^k^ik, := v "</• We interpret 



<[T3j)-(fT5j) as an algorithm iterating on the primal variables v and dual variables 71 over time 
t > as follows. For a generator bus j, given the current iterate (v(t),ir(t)), one can use 
'jjdL (v(t), n(t)) jdvj as the rate of change in Vj. For a load bus j, the current iterate i/j(t) is 
obtained as the solution of 

dL 



given n(t) and v~j(t) := (yi(t),i 7^ j). For a link G £ , given the current iterate (u(t), ir(t)), 
one can use —^ijdL (v(t), n(t)) /diXij as the rate of change in tt^-. 

In p3~])-(fl"5]), the stepsizes % and ^ can take any positive values, and the initial values 
(i/(0), 7r(0)) can also be taken arbitrarily. In particular, let 

7, = M-\ = 

and 

i/(0)=w(0), 7r(0) = P(0), 

where w(0) and P(0) are the frequency deviations and branch flow deviations at t = 0, the 
instant right after a step change P m in generation has occurred. We use continuous time t as 



the iterating time in the primal-dual algorithm (p~3|) — (f]~5|), and get a trajectory (u(t),ir(t)) for 
t > 0. We compare it with (u(t), P(t)), the trajectory of frequency deviations and branch flow 
deviations, and find that the the two trajectories take the same value for all t. Mathematically, 
(fT3|)-([T5|) is identical to ([5])-(|7]), if we identify v with u and rt with P, and take dj = dj(uj) in 
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For convenience, we collect here the system dynamics and load control: 

[dj + dj - P™ + Ff* - Pf\ for j eg (16) 

,m + pom _ p m for ■ e £ (17) 

f ) for(i,j)e£ (18) 

for j E J\f (19) 

for j G A/". (20) 





1 / 






= -Wj{ 









= Bij (ui 




= DjUj 


dj 


= k--v. 



The dynamics <[T6|)-([T9|) are automatically carried out by the power system while the active con- 
trol ( |20| ) needs to be implemented at each frequency-insensitive load. Let (d(t),d(t),uj(t), P(t)) 
denote a trajectory of frequency-insensitive loads, frequency-sensitive loads, frequency deviations 



and branch flow deviations, generated by the dynamics (16)-(20) of the load-controlled system 



Theorem 1. Every trajectory (d(t), d(t), ui(t), P(t)) generated by (|T6[) — (|20[) converges to a limit 
(d* , d* , cu* , P*) as t -> oo such that 

1) (d*,d*) is the unique vector of optimal load control for OLC; 

2) cu* is the unique vector of optimal frequency deviations for DOLC; 

3) P* is a vector of optimal branch flows for the dual of DOLC. 

We will prove Theorem [T] and other results in Section IV below. 



B. Implications 

Our main results have several important implications: 

1) Frequency-based load control: The frequency-insensitive loads can be controlled using 



their individual marginal cost functions according to ( |20| ), based only on frequency devi- 
ations tUj(t) (from their nominal values) that are measured at their local buses. Note that 
both the load control here and the generator droop control Q respond to the difference 
between the nominal frequency and the actual frequency, but they are complementary. Load 
control is activated immediately after a sudden generation-load imbalance because many 
loads can respond quickly. Droop control is slower than load control due to larger time 
constants associated with valves and prime movers, but it compensates for a large amount 
of generation-load imbalance to prevent the frequency from wandering outside the desired 
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limit. Optimal load control explicitly uses cost functions of individual loads to share the 
required load reduction/increase optimally among all loads in the network, whereas in 
droop control generator outputs are adjusted proportionally to frequency deviations so that 
the sum of implicit cost functions which are quadratic in the changes of generator outputs 
is minimized. 

2) Complete decentralization. The common operating frequency is a global signal that mea- 
sures the power imbalance across the entire network. Our result implies that the local 
frequency deviation 0Jj{t) at each bus turns out to convey exactly the right information 
about the global power imbalance for the loads themselves to make optimal decisions 
based on their own marginal cost functions. That is, with the right information, their local 
decisions turn out to be globally optimal. This result allows a completely decentralized 
solution without explicit communication among the buses. 



3) Reverse engineering of swing dynamics. The frequency-based load control ( |20| ) coupled 
with the dynamics (fT6|) — <fT9|) of swing equations and branch power flows serve as a 
distributed primal-dual algorithm to solve OLC and its dual DOLC. 

4) Frequency and branch flows. In the context of optimal load control, the frequency devi- 
ations Uj(t) emerge as the Lagrange multipliers of OLC that measure the cost of power 
imbalance, whereas the branch flow deviations Pij(t) emerge as the Lagrange multipliers 
of DOLC that measure the cost of frequency asynchronism. 

5) Uniqueness of solution. Lemma [2] implies that the optimal frequency to* is unique and 
hence the optimal load control (d*,d*) is unique. As we show below, the optimal branch 
flows P* are unique if and only if the network is a tree. Theorem [T] says nonetheless, 
that, even for mesh networks, any trajectory generated by the load control and system 
dynamics indeed converges to an optimal point, with the optimal value of P* dependent 
on the initial condition right after a change in generation. 

6) Optimal frequency. The structure of DOLC indicates that the frequencies at all the buses are 
synchronized at optimality even though they can be different during transient. However, the 
common frequency deviation cu* at optimality is in general nonzero. This implies that while 
frequency-based load control and the swing dynamics can resynchronize bus frequencies to 
a unique common value after a change in generation, the new frequency may be different 
from the common operating frequency before the change. To respect the tight frequency 
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regulation limits in power systems, the new steady-state frequency deviation should be 
made small. Simulations in Section [V] show that the new steady-state frequency deviations 
are reasonably small with OLC. Other mechanisms, such as isochronous generators 03 or 
automatic generation control [2], will be needed to drive the new operating frequency to 
its nominal value, through, e.g., integral control of the frequency deviation. 
Of course, many of these insights are well known; our results merely provide a fresh and unified 
interpretation within an optimization framework for frequency-based load control. 

IV. Convergence analysis 

This section is devoted to the proof of Theorem [T] and other properties as given by Theorems 
[2] - [3] below. Before going into the details, we first sketch out the key steps in establishing 
Theorem [1} the convergence of the trajectories generated by (fT6])-(|20]). 



1) Theorem [2j The set of optimal points (u*,P*) of DOLC and its dual and the set of 



equilibrium points of (fT6|)-(|20|) are nonempty and the same. Denote both of them by Z*. 

2) Theorem [3j If (A/", £) is a tree network, Z* is a singleton with a unique equilibrium point 
(co*,P*), otherwise (if (Af,£) is a mesh network), Z* has an uncountably infinite number 
(a subspace) of equilibria with the same to* but different P*. 

3) Theorem [TJ We use a Lyapunov-type technique to prove that every trajectory (u(t), P(t)) 



generated by (|T6|)-([20|) approaches a nonempty, compact subset Z + of Z* as t — > oo. 
Hence, if (A/", £) is a tree network, it is straightforward from Theorem [3] that any trajectory 
(w(t), P(t)) converges to the unique optimal point (to*, P*). If (Af, 8) is a mesh network, 
we show with a more careful argument that (u(t), P(t)) still converges to a point in Z + , 
as opposed to wandering around Z + . Theorem 1 then follows from Lemma [2j 
We now elaborate on these ideas. 



Given u, the optimal loads (d, d) are uniquely determined by (|T9])— ([20]) , hence we focus on the 
variables (co,P). Let C be the \J\f\ x\£\ incidence matrix with Cj e = 1 if e = (j, k) E 8 for some 
bus k E J\f, Cj e = — 1 if e = E 8 for some bus i E Af, and Cj e = otherwise. Recall that 
we assumed that the first \Q\ buses {1, . . . , \Q\} are generator buses and the remaining \C\ buses 
{\Q\ + 1, . . . , \Af\} are load buses. Decompose C into an \Q\ x \8\ submatrix Cg corresponding 
to generator buses and an \C\ x \£\ submatrix Cc corresponding to load buses, i.e., C = [ c ^]- 
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Similarly, let tug and tu c respectively denote the vector of frequency deviations at generator buses 
and load buses, so cu T = \ojg ux^\ . Let 

®0(uq) ■= and L g{ug,P) := $g(u g ) -UJ^CgP, 

jeg 

®c(u c ) := ^2$^) and L c (u c , P) := - Uc c cP- 



Identifying v with uj and n with P, we can rewrite the Lagrangian for DOLC defined in ( [12] ), 
in terms of uig and u>c, as 



L(U,P) = $M-W J CP = Lg(ug,P)+L C (u C ,P). 



Then (16)-(20) (equivalently, <pT3[>— ( 15 )) can be rewritten in the vector form as 



UJg 



dL c 



8Lg 



1 T 



dug 



due 



{ug) 



-C G P 



dur 



CrP 



8L 
dP 



u,P) 



~C T cu 



(21) 

(22) 

(23) 
(24) 



where Tg := diag^, j G Q) and E := diag(^-, E S). The differential algebraic equations 



(|22|)-(|24|) describe the dynamics of the power network when active load control is performed. 
A pair (u*, P*) is called a saddle point of L if 



L(u,P*) < L(co*,P*) < L(co*,P) for all 



(25) 



By [|T9l Section 5.4.2], (cj*,P*) is primal-dual optimal for DOLC and its dual if and only if 
it is a saddle point of L(co,P). The following theorem establishes the equivalence between the 



primal-dual optimal points and the equilibrium points of (22)-(24) 



Theorem 2. A point (u*,P*) is primal-dual optimal for DOLC and its dual if and only if it 
is an equilibrium point of (|22|)— ([24}) . Moreover, at least one primal-dual optimal point (u*,P*) 
exists and cu* is unique among all possible points (ui*,P*) that are primal-dual optimal. 

Proof: Recall that we identified v with u and n with P. In DOLC, the objective function $ 
is (strictly) concave over (by Lemma [I]), its constraints are linear, and a finite optimal u* 
is attained (by Lemma [2]). These facts imply that there is no duality gap between DOLC and its 
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dual, and there exists a dual optimal point P* |[T9l Section 5.2.3]. Moreover, (u*, P*) is optimal 
for DOLC and its dual if and only if the following Karush-Kuhn- Tucker (KKT) conditions [fT9l 
Section 5.5.3] are satisfied: 



Stationarity: 
Primal feasibility: 



LO; 



(CP*) T 
for all e £. 



(26) 
(27) 



On the other hand, (to*, P* 



(uig : uj£i P*) is an equilibrium point of (|22[) — (|24|) if and only if 

T 



dug 
diOr 



[OJ. 



T) 



CnP* 



[to, 



CrP* 



~C T u* = 0, 



which are identical to (|2"6])-(|2"7]). Hence, (u*,P*) is primal-dual optimal if and only if it is an 
equilibrium point of ([22])— ([24]) . The uniqueness of to* is given by Lemma [2j ■ 
From Lemma [2J we denote the unique optimal point of DOLC by = [^*\ 5 C ], where 

Iat G R m , \q E M) Gl and l c e M |£| have all their elements equal to 1. From ([26j)— ([27]) , define 



the nonempty set of equilibrium points of (|22|)-(|24|) (equivalently, primal-dual optimal points of 
DOLC and its dual) as 

"9$ 



(u,P) I u = CP 



(28) 



Let {u*ltf,P*) = (u*l g ,u*l c ,P*) G Z* be any equilibrium point of ((22j)-((24j). We consider a 
candidate Lyapunov function 



U(u,P) = -(ug-U*lg) T r g 1 (u J g 



E- % (P-P* 



(29) 



Obviously U(u, P) > for all (to, P) with equality if and only if ujg = u*lg and P = P*. We 
will show below that U(u,P) < for all (ui,P), where U denotes the derivative of U along 
the trajectory (u(t),P(t)). 



Even though U depends explicitly only on ujg and P, U depends on uc as well through ( |24| ). 
However, it will prove convenient to express U as a function of only tug and P. To this end, write 
< |23| ) as F(uj c ,P) = 0. Then ^(oj£,P) = ^f{wc) is nonsingular for all (uj c ,P) from the 
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proof of Lemma [T] in Section VIII- A By the inverse function theorem BUI , cu c can be written 
as a continuously differentiable function of P, denoted by lu c (P), with 



dP 



(P) 



du 2 c 



MP)) c c . 



-i 



Then we can rewrite the Lagrangian L(u,P) as a function of only (cog,P) as 
L(u,P) = Lg(ujg,P) + L c (u c (P),P) = L(ujg,P). 



We have the following lemma, proved in VIII-C, regarding the properties of L. 



Lemma 3. L is strictly concave in tug and convex in P. 



Rewrite (22)-(24) as 



OJg 



P = -E 



dL 

dbJg 

dL 

dP 



-1 T 



("G,P) 



Then the derivative of U along any trajectory (u(t),P(t)) generated by (|22|)— ([24]) is 



U(u,P) 



(O0g - U*lg) T Yg 1 O0g + (P- P*) T E^P 
— (Ug, P) (Cog - U*lg) - —(ojg, P) (P - P*) 

■L(ug,P) 



(30) 



(31) 



(32) 
(33) 



(34) 
(35) 
(36) 
(37) 



< L(ug,P)-L(u*lg,P)+L(ug,P* 
= L(uJg,CU*l c ,P*)-L(uJ*lg,P) 

< L(0J*l M ,P)-L(0J*lg,P) 
= Lg (C0*lg, P) + L C (CJ*1 C , P) - [Lg (iU*lg, P) + L C (cJ C (P), P)] 

< 0. (38) 
Here ((34]) follows from (f32|)— (|33]) . The inequality in p5] ) results from Lemma [3} The equality 



in d36J) holds since u c (P*) = uj*1 c by d26f . The inequality in (37 1 is due to L (u g , w*l c , P*) < 



L (oj*1m, P*) < L (oj*1m, P) by the saddle point condition ( f25| ). The inequality in ( [38] ) follows 
since uc(P) is the maximizer of L C (-,P) given P, by the concavity of L c in uj c and the 
definition of cuc(P)- 
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The next lemma, proved in Section |VIII-D[ characterizes the set in which the value of U does 
not change over time. 

Lemma 4. U(u, P) = if and only if either of the following two conditions holds. 
1) 



ug = uj*lg and CcP 



duj/ 



oo*h 



2) 



ug = u*lg and U1 C (P) = UJ*lc- 



Lemma [4] motivates the definition of the set 

E := {(u,P) \U(lu,P)=0} = \ (u, P) | u = u* \ N , C C P 



T 



(39) 



(40) 



(41) 



in which U — along any trajectory (u(t), P(t)). The definition of Z* in < [28| ) implies that 
Z* C E, as shown in Figure [2j As shown in the figure, E may contain points that are not in Z*. 




{a*0),P(0)) 



CO 



Fig. 2. E is the set on which U — 0, Z* is the set of equilibrium points of <|22[> — (|24[>, and Z + is a compact subset of Z* to 
which all solutions (u(t), P(t)) approach as t — > oo. Indeed, every solution (uj(t), P(t)) converges to a point (ui*,P*) G Z + 
that is dependent on the initial state. 



Nonetheless, every accumulation point (limit point of any convergent sequence sampled from 



the trajectory) of a solution (u(t), P(t)) of (|22[) — <|24]) is in Z*, as the next lemma shows 



Lemma 5. Every solution (u(t),P(t)) of ([22]) — (|24]> approaches a nonempty, compact subset 
(denoted Z + ) of Z* as t ->■ oo. 
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The proof of Lemma E\ is given in Section VIII-E The sets Z + C Z* C P are illustrated in 



Figure [2] Lemma [5] only guarantees that (u(t),P(t)) approaches Z + as t — > oo, but does not 
guarantee that it converges to any point in Z*. We now show that (u(t), P(t)) indeed converges 
to an equilibrium point in Z + . Indeed, the convergence is immediate in the special simple case 
when Z* is a singleton, but needs a more careful argument when Z* has multiple points. The 
next theorem reveals the relation between the number of points in Z* and the network topology. 

Theorem 3. 1) Suppose (M,£) is a tree, then Z* is a singleton. 

2) Suppose (M,£) is a mesh (i.e., contains a cycle if regarded as an undirected graph), then 
Z* has an uncountably infinite number of points with the same cu* but different P* . 



Proof: Recall that any point (u*,P*) G Z* is a solution of d26M27l>. Let h* := CP* = 
[|5 (to*)] T ■ Let (7 be the — 1) x \£\ reduced incidence matrix obtained from C by removing 
any one of its rows. Then C has a full row rank of \M\ — 1 lETI . Consider the corresponding 
equation 

CP* = h* (42) 

where h* is obtained from h* by removing the corresponding row. Since cu* is unique, so is h*. 
If (M,£) is a tree, then \£\ = \J\f \ — 1. Hence C is square and invertible, so P* is unique. If 
(N,£) is a (connected) mesh, then \£\ > \J\f\ — 1, so C has a nontrivial null space and there 
are uncountably many P* that solves ((42]). ■ 

With all the results above, we can now finish the proof of Theorem [TJ 

Proof of Theorem 1: For the case in which (N,£) is a tree, Lemma [5] and Theorem [3^1) 
guarantees that every trajectory(a;(t), P(t)) converges to the unique primal-dual optimal point 
(co*,P*) of DOLC and its dual, which, by Lemma [2j immediately implies Theorem [T] 

For the case in which (N,£) is a mesh, since U(cu,P) < for all (co,P), any solution 
(w(t), P(t)) for t > stays in the compact set {(cu, P)\U(u, P) < U(u(0), P(0))}. Hence there 
exists a convergent subsequence {(o;(tfe), P(tfc)), k E N}, where < t\ < t 2 < ••• and t k — > oo 
as k — > oo, such that lim^oo w(tfc) = co°° and limfc^oo P(tfe) = P°° for some (w 00 , P°°). Lemma 
IH implies that P°°) G Z + C Z*, and hence w 00 = w*Lv by ([28]). Recall that the Lyapunov 



function U in ( [29] ) can be defined in terms of any equilibrium point (u*, P*) G Z*. In particular, 
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select (cj*,P*) = (u/^P 00 ), i.e., 

U(u,P) := \{u g -un g ) T T g l {uj g -un g )+ l -{P-P™) T ~-\P-P™). 

Since U > and U < along all trajectories (u;(i), P(t)), U (u(t), P(t)) must converge as 
t — 7- oo. Moreover it converges to due to the continuity of U in both cu and P: 

lim U(u(t),P(t)) = lim U (u(t k ),P(t k )) = U{uf°,P°°) = 0. 

The equation above implies that the trajectory (u(t), P(t)) converges to (tu 00 ^ 00 ) E Z + C Z*, 
a primal-dual optimal point for DOLC and its dual. Theorem [T] then follows from Lemma [2j ■ 

Remark 2. The standard technique of using a Lyapunov function that is quadratic in primal-dual 
variables was first proposed by Arrow et al. [|22l . and has been revisited many times, e.g., in 
ll23ll [|24|. We apply a variation of this general technique to our particular problem and extend 
the results in the literature. First, with the algebraic equation ( [23] ) in the system, we take a 
Lyapunov function candidate that is quadratic in part of the primal variables cu g and the dual 
variables P, and show that it is indeed a Lyapunov function. Second, in the case when there are 
a sub space of equilibrium points due to the non-tree topology of the network, we show that the 
system trajectory converges to one of the equilibrium points instead of oscillating around the 
equilibrium set, without any modifications to the primal-dual algorithm like those in J24j. 

V. Case studies 

In this section, we illustrate the performance of OLC through simulation of the IEEE 68-bus 
New England/New York interconnection test system [3k The single line diagram of the 68-bus 
system is given in Figure [3} We run the simulation on Power System Toolbox lfT2ll . with two-axis 
subtransient reactance generator model, IEEE type DC1 exciter model, classical power system 
stabilizer (PSS) model, and AC power flow model on non-zero resistance lines. The detailed 
models and their parameters can be found in the data file and manual of the toolbox. 

In the system, there are 35 load buses serving different types of loads, including constant active 
current loads, constant impedance loads, and induction motor loads, with a total real power of 
18.23 GW. In addition, we add three loads to buses 1, 7 and 27, each making a step increase 
of real power by 1 pu (based on 100 MVA), as the change P™ in generation. We also select 
30 load buses, each having a load that is actively controlled based on OLC. In the simulation, 
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Fig. 3. Single line diagram of the 68-bus New England/New York test system. 



we take the same bounds [d, ~d\ with d = —d for each of the 30 controllable loads, and call 
the value of 30 x d the total size of controllable loads. We will show simulation results with 
different sizes of controllable loads below. In the simulation, the cost function of a controllable 
load dj is defined as Cj(dj) = d'j/(2a), with the same a over all the loads. As an example, 
we select a = 100 pu. To incorporate some practical consideration, the loads are not controlled 
continuously over time. Instead, they measure local frequencies and control their power every 
250 ms, which takes a relatively conservative estimate for the rate of load control 11251 . 

Since we have theoretically proved that OLC drives the system to a steady state where the cost 
of load control is minimized and total generation and total load are balanced, in the simulation we 
mainly focus on the transient performance, specifically, how the frequency and voltage change 
after a change in generation. We also look at how effective OLC is as a complement to the 
existing control mechanisms, such as the power system stabilizer (PSS), by enabling/disabling 



the PSS module in the simulation toolbox. Figures 4(a) and 4(b) respectively show the frequency 
and voltage at bus 66, under four cases: (i) no PSS, no OLC; (ii) with PSS, no OLC; (iii) no 
PSS, with OLC; and (iv) with PSS and OLC. In both cases (ii) and (iv), we take the total size 
of controllable loads as 1.5 pu. We observe that whether PSS is used or not, adding OLC can 
always improve the transient performance of frequency, in the sense that both the overshoot 
and the settling time (defined as the time after which the difference between the frequency and 
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Frequency at bus 66 



Voltage at bus 66 



no PSS. no OLC 

- - PSS, no OLC 
■--no PSS, OLC 
PSS, OLC 




0.992„ L 



no PSS, no OLC 
-PSS, no OLC 
-no PSS, OLC 
- PSS, OLC 




no OLC 



15 20 
Time (s) 



(a) 



(b) 



Fig. 4. The (a) frequency and (b) voltage at bus 66, under four cases: (i) no PSS, no OLC; (ii) with PSS, no OLC; (iii) no 
PSS, with OLC; (iv) with PSS and OLC. 



its new steady- state value never goes beyond 5% of the difference between its old and new 
steady-state values) are decreased. Using OLC also leads to a smaller steady-state frequency 
error. Comparison between cases (ii) and (iii) also suggests that using OLC solely without 
PSS produces a much better performance than using PSS solely without OLC. However the 
improvement of the transient performance of voltage is not as significant as frequency, which 
may be due to the fact that voltage depends more on reactive power injections while here OLC 
does not control the reactive power of loads. The effect of reactive load control to support voltage 
will be investigated in future work. 



To better quantify the performance improvement due to OLC we plot in Figures 5(a) - 5(c) 



the new steady-state value, the lowest value (which indicates overshoot) and the settling time of 
frequency at bus 66, against total size of controllable loads, as shown . Here PSS is always 
enabled. We observe that using OLC always leads to a higher new steady-state frequency 
(a smaller steady-state error), a higher lowest frequency (a smaller overshoot), and a shorter 
settling time, regardless of the total size of controllable loads. As the total size of controllable 
loads increases, the steady-state error and overshoot decrease almost linearly until a saturation 
around 1.5 pu. There is a similar trend for the settling time, though the linear dependence is 
only approximate. In summary, OLC improves both the steady-state value and the transient 
performance of frequency, and, in general, deploying more and larger controllable loads enables 
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Steady-state frequency at bus 66 




0.5 1 1.5 

Size of controllable loads (pu) 



it 59.96 



Lowest frequency at bus 66 



no OLC 
-OLC 



0.5 1 1.5 

Size of controllable loads (pu) 



(a) 



(b) 



Frequency settling time at bus 66 




0.5 1 1.5 

Size of controllable loads (pu) 



(c) 



Fig. 5. The (a) new steady-state value, (b) lowest value and (c) settling time of frequency at bus 66, against the total size of 
controllable loads. 



larger improvement. 

VI. Conclusion 

We have formulated an optimal load control (OLC) problem in power transmission networks 
where the objective is to minimize the cost of participation in load control subject to power 
balance across the network. We have shown that the dynamics of the swing equations and 
the branch power flows, coupled with a frequency-based load control, serve as a distributed 
primal-dual algorithm to solve the dual problem of OLC. Even though the system has multiple 
equilibrium points (nonunique branch power flows), we have proved that it nonetheless converges 
to an optimal point. Simulation of the IEEE 68-bus test system confirmed that the proposed 
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mechanism can resynchronize bus frequencies with significantly improved transient performance 
compared to using only local generator control mechanisms. 
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VII. Appendix: modeling details 
In this section, we first derive the branch flow dynamic model given by @-(|4]). Then, to show 



the accuracy of the model introduced in Section II-A (referred to as "the analytic model"), we 
use a more realistic simulation model developed in ll3l lfT2l (the same one we used in Section 
[V]) as a benchmark and demonstrate that the analytic model is a reasonable approximation of 
the simulation model. The key conclusions from simulations are summarized as follows: 

1) The internal and terminal voltage phase angles of the generator swing coherently, i.e., the 
rotating speed of a generator is always the same as the frequency at the generator bus. 

2) Different buses may have their own frequencies and buses that are far apart in electrical 
distance resynchronize at a similar timescale as the convergence time. 

3) The simulation model and the analytic model exhibit similar transient behaviors and steady 
state values of bus frequencies and branch power flows. 



A. Derivation of the branch flow model 

We assume that the system is always under the three-phase balanced condition, the frequency 
deviations Aujj are small, and the differences A9i — A9j between phase angle deviations are 
small across all the links G £. Specifically, Aujj is negligible compared to to , and an 

approximation of a quantity to the first order of A9i — A9j is reasonable. Now we show that the 
deviations APij in three-phase instantaneous power flows from their nominal values follow the 
dynamics in @-(|4]), by solving the differential equation that characterizes the line induction. 
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Without loss of generality suppose that buses i and j are wye-connected [|2) and each of the 
three lines has the same inductance L and zero resistance. Let the phase a voltages at buses i and 
j at time t be vf{t) = V2\Vi\ cos(w°t + 0? + A0i(t)) and ^ a (t) = v^l^l cos(a; t + 0j + A0 i (t)) 
respectively, and assume the voltage magnitudes are fixed. Denote the phase a current from % to 
j at time t by i%(t). 

For t < 0, suppose A0j(t) = for all the buses j. Hence the system is at a steady state with 
i%j(t) = a/2 I I Q | cos(o; t + 0°). From phasor calculations we have 

1/ I - N 
Mo I — 



0? = tan" 1 



1*5 


cos 0° - 




COS0,° 


1^ 


sin 0° - 




sin 0° 



where x i7 - := co°L, and |Vb| 



4-(0) = V2|J | cos 0; 



^l 2 + \Vj\ 2 ~ 2 1 Vi 1 cos(0° - 0°). Then, we have 
a/2 (I Vi I sin 0° - I Vj | sin 0j) 



X ; 



(43) 



For t > 0, we have 



whose solution is 



dt 



tW^W + l J o (vt(r)-v](r))dr 

«i£(0) + ^ [IVJI sin ( w °t + 0° + A0,(t)) - |y,| sin 0°] 
- [\ v i\ sin ("°* + #° + A W) - 1^1 sin ^°] 



(44) 



y/2 



X.; 



K| sin + 0° + A0i(t)) - I ^ I sin + 0° + A0,-(f))] 



where the approximate equality is due to the assumption that A0, = Auj are negligible compared 
to lu°, and the last equality is due to ( |43] ). 



May 6, 2013 



DRAFT 



25 



From (041) the instantaneous real power injection from i to j at phase a is 



Ptj 'T'",i 

Wi 



T-^- sin (u°t + 0° + A0 4 ) cos (w°t + 0? + A0i) 
- sin + 0° + A6j) cos + 0? + A0 4 ) 

1 >l ' |2 sin (2w t + 20? + 2A0 4 ) + sin (#o _ o + M _ A0 .) 



(45) 



mV] \ sin (2cA + 0° + 0° + A0, + A0,) 



Since we assumed the system is under the three-phase balanced condition, replacing 0° and 1 ? 



in (|45j) with 0° - |7T and 0° - |7T, we get replacing 0° and 0° in (|45]) with 0° + |7r and 
0° + |7r, we get p?j. Hence the three-phase instantaneous real power flow is (to the first order 

of A0i - A6j) 

WiWVA 



D ij ~ Pij Pij Pij 



py. + A p. . 



sin (0° - 0° + A0 ?; - A9j) 



(46) 



where 



is the nominal branch power flow, and 



A^ = 3^ 



cos(0°-0°)(A0 i - AOj) 



X{j 



is the deviation in branch power flow. By Aui = A0j, Awj = A0j, we get the branch flow 
dynamics in @-(|4]). 



B. Power flow behavior 

In Sections VII-B VH-C[ and VII-D we first simulate the IEEE 68-bus test system to a steady 
state (called "pre-change steady state") and then introduce the same step change in generation 
as in Section |Vj Then we compare the post-change behavior of the simulation with prediction 



of the analytic model introduced in Section II-A 
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We first check the branch flow dynamic model @-(|4]), which was derived above. Repeat it 
here as 

AP tj = B^Aui-Auj) for (i,j)eS (47) 

where Bij = 3\Vi\\Vj \ cos [9® — 9f) /xij is a constant under the assumption of constant voltage 
magnitudes and zero line resistances. We use the pre-change steady-state voltage magnitudes 
\Vi\, \Vj\ and angles 9®, 9® from the simulation to determine B^. Post change, we substitute 
the frequency deviations Auji(t) from the simulation into (j47]) to compute the trajectory Pij(t) 
that would result if APij is indeed proportional to the frequency difference, and compare it 
with the Pij(t) trajectory from the simulation. Since the simulation model is more detailed than 
the analytic model ([T])-((2]) for generator and load buses, using Auji(t) from the simulation to 
calculate Pij(t) isolates the behavior of branch flow modeled by ( |4"7| ). 

The results are shown in Figures 6(a) - |6(e) on five lines. The simulation model and the model 



given by < |47] ) exhibit similar transient behaviors and steady state values, suggesting that ( j47[ ) is 
a reasonable approximation of the simulation model. 

Specifically, to quantitative the accuracy of steady state, define the steady-state error as 

.j . . /modeled /simulated i nnO/ 

steady-state error := 1 — . x 100%, 

I /simulated | 

where "steady-state" refers to the post-change steady state, / mo deied refers to the steady-state 



frequency predicted by (j47J), and /simulated refers to the steady- state frequency given by the 
simulation. Figure [7] shows the errors on 84 of the 86 lines in the test system. Of these 84 
lines, all errors are within 11% and most within 2%. 

The two lines whose errors are not included in Figure [7J are the following: 

1) Line from buses 1 to 2: simulated power flow is -0.1 105 pu; modeled power flow is -0.0836 
pu; steady-state error is 24.35%. However the effect of the change in generation at steady 
state is similar between the simulation and the model prediction: the difference between 
pre-change steady state and the post-change steady state is 0.7258 pu in simulation and 
0.7527 pu from the model, yielding an error of 3.71%. 

2) Line from buses 52 to 42: simulated power flow is -0.0068 pu; modeled power flow is 
0.0198 pu; steady-state error is 393%. Note that the steady-state values themselves are 
much smaller than the average value (about 5pu). Moreover the difference between the 
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Real power flow on line 1 



Real power flow on line 3 
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Fig. 6. Real power flows on five of the lines in the 68-bus test system, both given by the simulation ("simulated" in the legend) 
and given by the model APij — Bij(AiUi — Acuj) ("modeled" in the legend). 
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Fig. 7. Histogram of steady-state real power flow errors (between the values given by simulation and predicted by the analytic 
model) on 84 out of the 86 lines in the IEEE 68-bus test system. 
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Fig. 8. Bus frequency is the same as the generator speed at bus 66 of the IEEE 68-bus test system. 

pre-change steady state and the post-change steady state is 0.950 pu in simulation and 
0.1215 pu from the model, yielding an error of 27.98%. 



C. Frequency behavior 



We check the assumption made in Section II-A that the internal and terminal voltage phase 
angles of the generator always differ by a constant, i.e., the rotating speed of a generator is 
always the same as the frequency at the generator bus. As an example, Figure [8] shows both the 
bus frequency and the generator speed at bus 66 in the IEEE 68-bus test system, which supports 
this assumption. 
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We then check a key modeling assumption that different buses may have their own frequencies 
and buses that are far apart in electrical distance resynchronize at a similar timescale as the 
convergence time. To this end we divide the 68 buses into the following 4 groups, with buses 
in each group being close in electrical distance to each other: 

1) Group 1 has buses 41, 42, 66, 67, 52, and 68; 

2) Group 2 has buses 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 
23, 24, 25, 26, 27, 28, 29, 53, 54, 55, 56, 57, 58, 59, 60, and 61; 

3) Group 3 has buses 1, 9, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 43, 44, 45, 46, 47, 48, 
49, 51, 62, 63, 64, and 65; 

4) Group 4 has bus 50 only. 

Figure [9] shows all 68 bus frequencies from the simulation, in four groups. We see that the 
frequencies at different buses within the same group are almost identical even during transient, 
but the bus frequencies of different groups are different during transient. Moreover the time it 
takes for these different frequencies to converge to a common system frequency is on the same 
order as the time for these frequencies to reach their (common) equilibrium value. 



D. Accuracy of the analytic model 

To show the accuracy of the analytic model <[l])-(|4]), we model the 68-bus test system as a 
4-node network, where each node represents an area (group of buses shown in Figure [9]) within 
which all buses have roughly the same frequency. Of these 4 nodes, we model nodes 1-3 as 
generator nodes with positive inertia by the swing dynamics ([T]), and model node 4 (which has 
bus 50 only) as a load node with zero inertia by the algebraic equation ([2]). 



Figure 10 shows the frequencies at four of the 68 buses, which are representatives for the 



four groups, both given by the simulation and given by the analytic model. Figure 1 1 shows the 
real power flows on two of the lines connecting two groups, both given by the simulation and 
given by the analytic model. As we can see, the analytic model is a reasonable approximation 
of the more realistic simulation model. 
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Fig. 9. Bus frequencies in each of the four groups. 



VIII. Appendix: proofs 



A. Proof of Lemma [7] 



From ( fTT| ), either c'-{dj{u)) = v or d'j(v) = 0. Hence, in ( [10] ) we have 



dv 



(c^djiy)) - udj(y)) = c^d^d'^v) - djiy) - vd'^v) = -dj(v), 



and thus 



<9$ 



[V] 



-djM-DM + Pp. 



Hence the Hessian of $ is diagonal. Moreover, since dj{yf) given by (11) is nondecreasing in 



za , we have 



-d'Avi 



Q V 2 VI 3 



Da < 0. 
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Fig. 10. Frequencies at four buses which are the representatives for the four groups, both given by the simulation ("PST 
simulated" in the legend) and given by the analytic model ("modeled" in the legend). 



Hence $ is strictly concave over M) \. □ 
B. Proof of Lemma [2] 

Let g denote the objective function of OLC with the domain V :— [t^, di\ x • • • x [dun, dpv[] x 
IRl^l. Since Cj is continuous on [dj, df\, ^ ■ Cj(dj) is lower bounded, i.e., J2j c j(dj) > C_ for some 
C > — oo. Let (d', d!) be a feasible point of OLC (which exists by Condition [T|). Define the set 
V ■= | (d, d\ e V d) < 2Dj(g(d', d') - C) for all j e Af\. Note that for any (d, d\ e V\V, 
there is some i E J\f such that df > 2Di(g(d',d') - C), thus 

9(dJ)>C+^>g(d',d'). 
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Hence, any optimal point of OLC must lie in V . By Condition [T] the objective function g of 
OLC is continuous and strictly convex over the compact, convex set V, thus it has a minimum 
g* > — oo attained at a unique point (d*, d*) G V . 

Note that OLC has a feasible point (d, d) G relint V, where relint V denotes the relative 
interior of V [fT9ll . Indeed, let (<f, d') G T> be a feasible point of OLC, then we can obtain (d, d) 
by letting dj = (c^ + dj) /2, dj = d'j — dj+d'j. Also note that the only constraint of OLC is affine. 
Hence, there is zero duality gap between OLC and its dual, and a dual optimal u* is attained since 



g* > -oo EE Section 5.2.3]. By Section VIII-A £ i6A ^'0) = - (d'j{v) + Dj) < 0, 



i.e., the objective function of the dual problem of OLC is strictly concave over R, which implies 
the uniqueness of v* . Moreover, by the definition of dual problem, the optimal point {d* ,d*^ 
of OLC should satisfy d* = dj{v*) given by dTTb and d* = Djis* for all j G M. □ 



C. Proof of Lemma [3] 

From the proof of Lemma |IJ the Hessian ^-(ug,P) = ^jr-i^g) is diagonal and negative 



definite for all ug G Therefore L is strictly concave in tug. Moreover, from pTj ) and the 
fact that §gr {u c {P), P) = 0, we have 

BL 

— (Ug, P) = -U T gCg - U T C (P)C C . (48) 
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Therefore we have (using ( |30| )) 

2 L 



r T duJ C 

~ Cc ~dP 



(P) 



-cl 



MP)) 



-1 



Cr, 



9 <4 



is diagonal and negative definite. Hence jp%(ujg, P) is positive 



From the proof of Lemma jlj 
semidefinite and L is convex in P (L may not be strictly convex in P because Cc is not 
necessarily of full rank). □ 

D. Proof of Lemma [?] 

The equivalence of ((40]) and ((39]) follows directly from the definition of uc(P). To prove that 



((40]) is necessary and sufficient for U(u, P) = 0, we first claim that the discussion preceding 
the lemma implies that (u, P) = (u)g,ujc,P) satisfies U(u,P) = if and only if 

dl 



cog 



,fl g and — (u g ,P)(P-P*) 



0. 



(49) 



Indeed, if (49) holds, then the expression in (34) evaluates to zero. Conversely, if U(u, P) = 0, 



then the inequality in (35) must hold with equality, which is possible only if tug = co*lg since L 
is strictly concave in tug. Then we must have ^(cug, P) (P — P*) =0 since the expression in 



( p4| ) needs to be zero. Hence we only need to establish the equivalence of ((49]) and ( |40| ). Indeed, 

with tug = u*lg, the other part of ( |49] ) becomes 

dl 



dP 



yi g ,p) (p-p*: 



[u,% u T c (p)]c(p-n 

[0 ujI(P)-u*1 t c ]C(P-P*) 
(u c (P)-u*l c ) T C c (P-P*) 



(u c (P)-u*l c ) 



r 



MP)) 



{un c ) 



(50) 
(51) 

(52) 



due Oujc 

where ((50) follows from @, (j5T) follows from 1^0 = 0, and ([52) follows from ((23]) and @. 
Note that $£ is separable over w,- for j e £ and, from ( [10] ), 3>i-(a>j) = —dj(ujj 
Writing £>£ := diag (/}.,, j G £), we hence have 



pm 
3 ' 



dL_ 
dP 



U*lg,P) (P-P*) 



(co c (P) -LU*l c ) T D c (lu c (P)-uj*Ic) 



(53) 



+ MP) - MP)) - djM) ■ 

Since dj(tOj) given by ( fTT| ) is nondecreasing in cuj, each term in the summation above is 



nonnegative for all P. Hence ((53]) evaluates to zero if and only if loc(P) = ui*lc, establishing 
the equivalence of ((49]) and ((40]). □ 
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E. Proof of Lemma [3] 

The proof of LaSalle's invariance principle in 11261 Theorem 3.4] shows that (u(t),P(t)) 
approaches its positive limit set Z + which is nonempty, compact, invariant, and a subset of 
E, as t — > oo. It is sufficient to show that Z + C Z* , i.e., considering any point (ui,P) = 



(ug,uc,P) E Z + , to show that (u,P) E Z*. By fl28| ), (|4TJ) and the fact that (w,P) G £, we 
only need to show that 

~d<f>g 



C C P 



1 T 



due 



(ug) 



(54) 



Since Z + is invariant with respect to ([22])— ([24]), a trajectory (u(t),P(t)) that starts in Z + must 



stay in Z + , and hence stay in E. By ( |4T| ), wg(t) = for all £ > 0, and therefore u)g(t) = 



for all t > 0. Hence, by (22), any trajectory (uj(t), P(t)) in Z + should satisfy 



C P(t) 



9$c 



T 



for all t > 0, which implies that (J54l) holds for any (u,P) E Z 



□ 
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